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Introduction 

The ability of near-infrared (NIR) light-based techniques to noninvasively image 
and analyze tissue structure and function promises their great potential for detection 
and diagnosis of breast cancer. Optical diagnostic techniques allow us to not only 
enhance the existing capabilities, but to eliminate the need for physical biopsies. In 
addition, optical imaging is inexpensive and portable, which indicates that optical 
imaging could be an ideal candidate for routine breast screening. However, since the 
scattering properties of tissues convolute re-emitted NIR signals, the extraction of 
pertinent information continues to remain elusive. An understanding of light propagation 
and light-tissue interaction is required before the optical technologies can substantially 
impact diagnostic medicine. 

Research efforts in the Biomedical Optics Laboratory at Clemson University are 
focused on the biophysics of light propagation and light-tissue interaction in order to 
engineer appropriate approaches for noninvasive breast imaging and spectroscopy. 
Specifically, we are developing indirect optical/fluorescence approaches using photon 
migration measurements in the continuous-wave and frequency domains. These 
indirect optical/fluorescence approaches or image reconstructions are computationally 
based on the powerful finite element methods. A CCD-based optical spectroscopic 
imaging system is already operational in our laboratory for continuous-wave 
tomographic photon migration measurements, while a frequency-domain system is still 
under construction. Using these optical systems coupled with our finite element based 
reconstruction algorithms, we will be able to extract spatial/spectroscopic maps of tissue 
optical properties, lifetime and/or yield of endogenous and exogenous fluorescent 
probes. Since metabolic tissue states can be identified by our approaches, diagnostic 
information is also obtained in addition to detection of tumor. 

This Career Development application for support of Dr. Huabei Jiang will facilitate 
the establishment/continuation of these research activities. Interdisciplinary interactions 
with the Greenville Hospital System (Greenville, SC) will be enhanced, which insures 
the direction of research towards a clinically pertinent and feasible system. 

Body 


This report describes work accomplished during the first year of a proposed four- 
year study. This Career Award supports Dr. Jiang’s research on optical and 
fluorescence imaging using both continuous-wave and frequency-domain 
measurements. The focus of the proposed work in Year 1 is the implementation of 
some enhancing schemes in our existing 2D reconstruction codes as well as some 
preparation for phantom studies. This will pave the way for us to conduct 2D image 
reconstruction in Year 2 for simple experimental configurations with dye-free phantom 
background. 

Software work: while we have implemented some of the proposed image enhancing 
schemes in the current 2D codes including the total variation minimization, weighted 
least squares criterion and low pass filtering, we have also developed two additional 
novel schemes that were not proposed originally. Our existing algorithms require an 
additional calibration measurement with a homogeneous phantom in order to determine 


1 



H. Jiang/DOD Annual Report 

the excitation source strength and the boundary conditions coefficient that are critical for 
a successful reconstruction. The first scheme, which uses the idea of normalizing the 
photon density in the reconstruction algorithm, allows for the reconstruction of optical 
property images without measuring the excitation source strength. The second scheme, 
which is based on a simple least-squares minimization between the measured and 
computed photon densities at the boundary, can provide us the boundary conditions 
coefficient. The normalizing scheme-based algorithm eliminates the need of absolute 
measurement data for reconstruction, yet provides us absolute quantitative optical 
images. We have tested these two schemes using simulations in which noise-free and 
noisy “measured” data are applied. We have also confirmed these simulations using 
tissue-like phantom experiments (see Phantom Experiments section below). For the 
detail of the two enhancing schemes and the simulation and experimental results, see a 
manuscript entitled “Enhancing Diffusion Optical Tomography through A Normalizing 
Reconstruction Algorithm and A Simple Least-Squares Minimization Scheme” provided 
in the Appendix to the Summary Report. 

Hardware Development: in order to conduct the proposed fluorescence lifetime 
imaging work, we must first construct a frequency-domain optical system. Support from 
NIH and Greenville Hospital System has allowed us to complete this system. We have 
built a 16x16 channel automatic frequency-domain imaging system that allows us to 
collect tomographic data within eight minutes at the present time. The system has been 
described in detail in a second manuscript entitled “Development of a combined optical 
and fluorescence imaging system in frequency-domain for breast cancer detection” 
provided in the Appendix to the Summary Report. 

Phantom Experiments: using the developed frequency-domain imaging system, we 
have conducted a series of phantom experiments to confirm the two new image 
enhancing schemes. The successful results have been presented in the first manuscript 
included in the Appendix. In addition, we have begun to conduct phantom experiments 
for fluorescence lifetime reconstruction which was proposed to take place during the 
second year of the project. In these experiments, we present successful quantitative 
simultaneous reconstruction of both lifetime and yield images from frequency-domain 
measurements using indocyanine green (ICG) dye and tissue-like phantoms. To our 
knowledge, the results shown here are the first experimentally reconstructed images for 
which simultaneous recovery of both lifetime and yield profiles has been achieved with 
absolute ac excitation and fluorescence data. We show the ability of our algorithm to 
reconstruct images in which single- and multi-target are embedded in a circular 
background. Our experimental setup used was the automated multi-channel frequency- 
domain system mentioned above. The system employed a radio-frequency intensity- 
modulated near-infrared beam. The laser beam at 785 nm was sent to the phantom by 
16 fiber optic bundles coupled with a high precision moving stage. The diffused 
radiation was received by another 16 channel fiber optic bundles and delivered to a 
PMT. A second PMT was used to record the reference signal. These PMTs were 
supplied at a radio-frequency modulated current with 0.1-1 KHz shift. The intermediary 
frequency signal obtained from the PMTs was processed using a National Instruments 
board. To increase the dynamic range of intensities that the PMT can detect, an 
automated filter wheel with pre-calibrated neutral density filters was added to the 
system. For every source position, 16 measurements for each detector were made, 
taking alternatively 100 ms samples for sample and reference signals. Fluorescence 
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signals were obtained through an 830nm interference filter placed in front of the 
detection PMT. ac intensity and phase shift between reference and sample signals were 
obtained using FFT Labview routines. The total data collection time for 256 
measurements was 8 minutes. 

In our experiments we used a 50 mm diameter solid phantom (0.5% Intralipid 
+India ink+Agar) as the background medium, making the absorption coefficient 
|i a =0.0005/mm, and the reduced scattering coefficient |i' =0.5/mm. Both single- and 
two-target configurations were used. The solid targets contained ICG dye with different 
concentrations (0.5 or 1 |iM) plus 0.5% Intralipid and Agar (perfect uptake of the dye 
was considered here, i.e., no dye was present in the background). The single target had 
12.6 or 7.8 mm in diameter, while both targets had 10.3 mm for the two-target case. The 
2D finite-element mesh used had 249 nodes and 448 elements for both forward and 
inverse solutions. The image results have been presented in a third manuscript entitled 
“Fluorescence lifetime tomography of heterogeneous turbid media using frequency- 
domain measurements” provided in the Appendix to this Summary Report. 

We have shown that we are able to obtain simultaneous reconstruction of both 
lifetime and yield images in turbid media using ac excitation and fluorescent data. While 
an ideal, perfect uptake of fluorescent dye in the target has been assumed, this study 
has clearly demonstrated the feasibility of fluorescence lifetime imaging in turbid media 
using the reconstruction approach described here. To add fluorophores into the 
background would primarily challenge our hardware system since the fluorescent 
signals would be weakened in this case. However, we do not anticipate considerable 
difficulty to obtain enough signal-to-noise ratio for reconstruction when the fluorophores 
are present in the background since the dye uptake ratio between the tumor and normal 
tissue can be as high as 10:1. In this regard, we plan to conduct experiments in Years 
2-3 of this project. 

Key Research Accomplishments 

1. We have developed a number of novel schemes that can enhance our current 2D 
reconstruction software. 

2. We have constructed and tested a muti-channel frequency-domain imaging system. 

3. We have conducted phantom experiments that confirmed our software 
enhancement. We have also performed successful fluorescence phantom studies 
that were proposed to occur during the second year of the project. 

Reportable Outcomes (see the Appendix to this Summary Report) 

Conclusions 

We have made a significant progress that has exceeded the statement of work 
proposed for Year 1 of this project. Given the successful first year, we will be able to 
fulfil or exceed the work statement for Year 2, and we will carry these successful works 
over into the final year of this project. 
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Abstract 


In this paper, we present two schemes that can enhance our established finite-element 
based optical image reconstruction algorithms. Our existing algorithms require an additional 
calibration measurement with a homogeneous phantom in order to determine the excitation 
source strength and the boundary conditions coefficient that are critical for a successful 
reconstruction. The first scheme, which uses the idea of normalizing the photon density in the 
reconstruction algorithm, allows for the reconstruction of optical property images without 
measuring the excitation source strength. The second scheme, which is based on a simple least- 
squares minimization between the measured and computed photon densities at the boundary, can 
provide us the boundary conditions coefficient. The normalizing scheme-based algorithm 
eliminates the need of absolute measurement data for reconstruction, yet provides us absolute 
quantitative optical images. We test these two schemes using simulations in which noise-free and 
noisy “measured” data are applied. We further confirm these simulations using tissue-like 
phantom experiments. 
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1. Introduction 


There has been considerable interest in the application of diffusion optical tomography 
for tissue diagnostics as evidenced by recent feature issues on this topic 1 ' 3 . Among various 
imaging approaches under developing, diffusion optical tomography based on reconstruction 
algorithms has received particular attention since it can offer both structural and functional 
information of tissue 4 ' 17 . This unique capability may allow for non-invasively differentiating 
between benign and malignant tumor, for example. In this regard, the core task is to develop an 
effective reconstruction algorithm for forming an image of the spatial distribution of tissue 
optical absorption and scattering properties. 

We have presented a finite-element based reconstruction algorithm 18 , which has been 
evaluated using extensive simulated and experimental data in the last few years 7,18 ' 22 . These 
studies clearly suggest that we have identified a powerful framework with which to pursue 
quantitative image reconstruction of both absorption and scattering property profiles. A 
calibration procedure was used in all of the experimental imaging studies performed, however. In 
this procedure, we used homogeneous medium measurements to calibrate the source term and 
the boundary conditions (BCs) 7,19 . Although it was straightforward to do it, a trial-and-error 
process was needed to complete this calibration procedure. We found that the amplitude of the 
source largely determined the computed photon density amplitude and BCs could affect the 
computed overall light distribution considerably, which suggests that this time-consuming 
procedure may introduce significant errors to the reconstruction if not carefully done. In 
particular, this calibration procedure may not be viable for in vivo data 22 . 

In this paper, we present two schemes that can enhance our established finite-element 
based optical image reconstruction algorithms. The first scheme, which uses the idea of 
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normalizing the photon density in the reconstruction algorithm, allows for the reconstruction of 
optical property images without measuring the excitation source strength. The second scheme, 
which is based on a simple least-squares minimization between the measured and computed 
photon densities at the boundary, can provide us the boundary conditions coefficient. The 
normalizing scheme-based algorithm eliminates the need of absolute measurement data for 
reconstruction, yet provides us absolute quantitative optical images. We perform numerical 
simulations to test the two schemes where noise-free and noisy “measured” data are applied. We 
confirm these simulations using tissue-like phantom experiments. 

2. Normalizing Scheme 

Our image reconstruction algorithm has been based on the 2D diffusion 
approximation 18,19 . For cw light illumination, the steady-state diffusion equation can be stated as 
V • D(r)V<D(r) - |i.(r)0(r) = ~S(r) (1) 

where <b(r) is the photon density, D(r) is the diffusion coefficient, and |l a (r) is the absorption 
coefficient. The diffusion coefficient can be written as 


D(r) = 


1 

3[fX a (r) -t- (J.; (r)] 


( 2 ) 


where (l' s (r) is the reduced scattering coefficient, S 0 is the source term in (1) which for a point 
source can be written as S = S 0 S(r-r 0 ), where S 0 is the source strength and 8(r-r 0 ) is the 
Dirac delta function for a source at r 0 . 

2.1. Non-normalizing Reconstruction Algorithm 


Our non-normalizing reconstruction algorithm (in short, Non-NRA) has been described 
in details in Ref. 18. Here we just provide an overview in order to present the normalizing 





reconstruction algorithm (in short, NRA). Following the procedures outlined by Paulsen and 
Jiang 18 , we can obtain a finite element discretization of (1) 


(A-aB)&=C 


where the elements of matrix A are a^ = /-^D k <|) k V<|)j-V^-^M^<)>£<t>j<t>; ) where ( ) 


r 

indicates integration over the problem domain; the elements of matrix B are by = l^j^ds 

j=i 

where <E> = [<3> 1 ,^> 2 ,--*-<I> N ] T , is the photon density at node i, N is the number of nodes in 
the finite elements mesh, M is the number of boundary measurements, j> expresses integration 


over the boundary surface where type m boundary conditions (BCs), -DVd> • n = a<E> , have 
been applied; a is a parameter relating to the internal reflection at the boundary; the elements of 
vector C are In (3), O, D and (l a have been expanded as the sum of coefficients 

multiplied by a set of locally spatially varying Lagrange polynomial basis functions <j)j ,<t>^ • 

Note that the expansions used to represent the diffusion and absorption coefficient profiles in (3) 
are K and L terms long where K * L * N in general; however, in the work reported here K = L 
= N. 

Based on a Newton's iterative scheme as described by Paulsen and Jiang , the following 
matrix equation can be obtained for providing the inverse solution of interest 

(3 t 3 + Xl^x = 3 T [<F (m) -<D (C) ] (4) 

where Ax = (AD a ,AD 2 ,->",AD N Ap a 1 ,Ap. aj2 ,----,A|J, a N ) T is the updates of the diffusion and 

absorption coefficient profiles; x represents either D or (i a ; 3 is the Jacobian matrix which 

consists of derivatives of the photon density with respect to % and 3 T is the transpose of 3 , 
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0^"9 and 0^ are measured and computed photon densities respectively; X is a regularization 

parameter which is used to regularize or stabilize the decomposition of 3 X 3 and I is a unit 
matrix. 


The elements of the Jacobian matrix are computed from the following matrix 
relationship, which is the differentiation of equation (3) 


(A-aB) 


90 = 9c 
9 % 9 % 


9A 

9% 


(5) 


Evaluation of (4) for the update vector % requires the solutions of (3) and (5), which 
depend on the current estimate of %. Once the update vector has been computed, a new % vector 
becomes available and the procedure is repeated until a converged value of the difference norm 
between the measured and computed data is reached. Note that, when reconstructing 2D images 
from experimental data using this Non-NRA, both the source term, S, and the boundary 
conditions coefficient, a, need to be determined through a calibration procedure with a 
homogenous phantom 7,19 . 

2.2. Normalizing Algorithm 

Based on (3) and (5), we derive a normalizing reconstruction algorithm. Assume that 0 O 
is a known measurement at a node on the boundary. Dividing 0 O at the both sides of equation 
(3), we can obtain a normalized form of equation (3) 


(A-aB)—= — 

$0 


Similarly, a normalized differentiation of matrix equation (5) leads to 

9 


(A-aB) 


9x 


(4>'| 

1 

f - \ 

9A 1 9C _ 9 
— 1 c 

f 1 1 


T 

oT 
^ °) 

9% 0 O 9% 9% 



( 6 ) 


(7) 
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(7) can be rearranged as follows 



By comparing (3) and (5) to (6) and (8), we can easily obtain the elements of the 
normalized Jacobian’s matrix which are needed to construct the normalized form of (4). Since 
the equations (5) and (8) originate from the same diffusion equation (3), both non-normalizing 
and normalizing algorithms should provide the same reconstruction of the optical properties. The 
main advantage of NRA over Non-NRA is that the former is independent of source intensity. In 
other words, we need not calibrate the source term in measurements under experimental or 
clinical conditions, which usually is a time-consuming procedure and may cause additional 
errors in measured data. In Section 4, we will confirm our NRA scheme by using simulations and 
experiments. 

3. Least-Squares Minimization Scheme 

We use a simple least-squares minimization scheme to determine the boundary 
conditions coefficient, a. That is, we compute the following X 2 value as a function of a, 

M 

X 2 =£(<D[ m) -^ c) ) 2 (9) 

i=l 

where M is the number of boundary measurements; is the measured photon density from a 
given experimental inhomogeneous medium (a phantom or in vivo tissue), and 3>[ c) is the 
computed photon density from a numerical simulation homogeneous medium with the same 
geometry as the experimental medium. We argue that the minimum of X 2 corresponds to the 
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correct value of a associated with the experimental inhomogeneous medium. This argument is 
confirmed in the following section using simulations and experiments. 

4. Experimental Materials and Methods 

Our experimental setup is an automated multi-channel frequency-domain system that has 
been described in detail elsewhere. 23 In this system, a radio-frequency intensity-modulated near- 
infrared beam at 785 nm is sent to the phantom by 16 fiber optic bundles coupled with a high 
precision moving stage. The diffused light is received by another 16 channel fiber optic bundles 
and delivered to a PMT. A second PMT is used to record the reference signal. These PMTs are 
supplied at a radio-frequency modulated current with 0.1-lKHz shift. The intermediary 
frequency signal obtained from the PMTs is processed using a National Instruments board. To 
increase the dynamic range of intensities that the PMT can detect, an automated filter wheel with 
pre-calibrated neutral density filters is added to the system. For every source position, 16 
measurements for each detector are made, taking alternatively 100 ms samples for sample and 
reference signals, dc, ac and phase shift between reference and sample signals are obtained using 
FFT Labview routines. The total data collection time for 256 measurements is 8 minutes. We just 
used measured dc data to reconstruct the absorption and scattering images in this study. 

In our experiments we used a 50 mm diameter cylindrical solid phantom (1% 
Intralipid+India ink+Agar) as the background medium, making the absorption coefficient, 
fi a =0.005/mm and the reduced scattering coefficient, fi' =1.0/mm. Both single- and two- 

target configurations were used. In the single-target case, a target with 15.0 mm in diameter was 
examined. In the two-target case, one 8.0 mm and one 12.5 mm diameter target were embedded 
in the background. |i a of the target varied from 0.020/mm to 0.040/mm; |i' of the target was the 
same as that of the background. These values of the background and target optical properties are 
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similar to those for normal and tumor breast tissues. 24 The 2D finite-element mesh used had 249 
nodes and 448 elements for both forward and inverse solutions. Note that here we do not need to 
use the dual mesh scheme as used before 21,22 in order to obtain quantitative reconstruction by dc 
data. All the reconstructions were performed in a 400 MHz Pentium II PC. 

5. Results 

In this section, we will evaluate the above NRA and the simple least-squares 
minimization scheme using simulated and experimental data. The parameters used in the 
experiments have been given in Section 4. For simulations, the background geometry and optical 
properties were exactly the same as that for experiments. A single target with a radius of 10 mm 
was embedded in the background. The optical properties in the target were ji a =0.01/mm and 

|i’ = 2.0/mm. The mesh used had 241 nodes and 460 triangular elements. A circular symmetric 

configuration of 16 source and 16 detectors were used, providing a set of 256 "measured" data. 
Assume that the data is noise free and has 5% noise respectively, we have tested both NRA and 
Non-NRA for comparison. For the reconstructions using NRA, we selected the maximum value 
of the photon densities at the boundary to normalize each set of "measured" data. To determine 
the boundary conditions coefficient, a, we computed the X 2 as function of a. 

Figs. 1 and 2 show the reconstructed images using Non-NRA and NRA in which 0% and 
5% noise were added to the “measured” data, respectively. When different source intensity 
values were used, all the images for both p. a and fi' reconstructed using NRA scheme are 

exactly the same as the images shown in Fig. 1(c) and (d) under the condition of 0% noise and in 
Fig. 2 (c) and (d) under the condition of 5% noise. For this reason, these images with different 
source intensities are not presented here. Table 1 gives the quantitative information about the 
RMS errors for the images shown in Figs. 1-2. The image RMS errors are defined as 
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V 1/N Si=,t(5Cexa« — X reconstructed ) ^ X exact P where N is the number of sample positions and x is 
either (i a or ji'. 

In Fig. 3 (a), X 2 is depicted as a function of a without noise added to the “measured” 
data when the optical properties used for computing <3>[ c) are the same as the background values, 
2%-off the background values, and 20%-off the background values, respectively. Fig. 3 (b) 
presents the same results as Fig. 3 (a) except that 5% noise was added to the “measured” data. 
Note that the exact value of a used here is 0.7. 

Fig. 4 displays X 2 versus a for the two sets of experimental data. As can be seen, a 
minimum value of the X 2 error is clearly reached at a = 0.57 for the single target case and at 
a = 0.52 for the two target case. Using these derived values of a, the absorption and scattering 
images were reconstructed and are shown in Fig. 5 for both experimental configurations. As 
evidenced, the images are clearly recovered. 

6. Discussion and Conclusions 

The results presented in the previous section confirm the NRA scheme outlined in 
Section 2. It is clear that the reconstructions of both (i a and |i' images using NRA are 
independent of the excitation source intensity. As can be seen from Figs.l and 2, when no noise 
was added to the "measured" data, the visual quality, in terms of the target shape and size, of (i a 

images reconstructed using Non-NRA and NRA is about the same, while the visual quality of p.' 

image recovered by NRA is even better than that by Non-NRA. Most interestingly, when 5% 
noise was added to the "measured" data, we obtained much better reconstructed images for both 
p. a and p,' using NRA than that using Non-NRA in terms of the target shape and size recovered. 
This suggests that NRA is less sensitive to the noise effect than Non-NRA for reconstructing 
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both jx a and p.' images. It is also interesting to note that the overall recovery of the target shape 
and size for ]i a images is better than that for |i' images, no matter Non-NRA or NRA is used. 

Table 1 presents the average values and the RMS errors of optical properties recovered. 
In terms of the average optical property values, Non-NRA and NRA gave similar accuracies in 
most cases except that NRA provided a better value of p,' than Non-NRA for the target under 
the condition of 5% noise (p,'=1.64 from NRA whereas p,'=1.27 from Non-NRA). In terms of 
the RMS errors, overall Non-NRA provided slightly better results than NRA for both p, a and p' 

under the condition of 0% noise. However, under the condition of 5% noise, neither Non-NRA 
nor NRA showed consistent better results than its counterpart. The relative high RMS errors in 
some of the cases here were caused by the relative high or low values of (l a or p,' at some 

locations. These relatively high or low values of |i a or p a can be easily smoothed by using a low 
pass filtering 19 ’ 21 , which would provide much lower RMS errors. 

No matter 0% or 5% was added to the “measured” data, Fig. 3 shows that X 2 reaches its 
single minimum right at the exact value of a =0.7, when the optical property values used for 
computing d>[ c) are the same as the background, or just off 2%. The minimum occurs at a =0.64 

when the optical property values used for computing 0[ c) are off 20% relative to that of the 
background under both the conditions of 0% and 5% noise. Clearly the amount of noise added to 
the “measured” data does not appear to affect the determination of a. The optical property 

values used for computing 0[ c) are the single important factor impacting this determination. 
Given the fact that the optical properties for the normal background tissue as a priori are 
typically known with an accuracy better than 10% error 24 ' 26 , we argue that the method described 
here offers a viable way for determining a in an in vivo setting. We have found that a significant 
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shift of a value (about a =0.45) presented when the optical property values used for computing 
0[ c) are 100%-off the background optical property values. 

From experimental data. Fig. 4 shows that X 2 reaches the single minimum at a =0.57 
and a =0.52 for the single and two target cases, respectively. The exact values of a are not 
available since it is difficult for us to obtain such values given our experimental conditions. 
Theoretically these values would depend upon the ratio of the refractive indexes of the phantom 
or tissue and air or surrounding medium. 27,28 In our experiments, most of the boundary space was 
surrounded by the thirty-two 3 mm source/detection optic fibers. This means that the surrounding 
medium used in our study was a combination of the optic fibers and air (about 75% of the 
medium was fibers and 25% was air). Since the refractive index of the glass optic fibers is about 
1.53, this would give us an effective index of 1.40 for the mixed medium which is larger than 
that of the phantom used (1.33 for the Intralipid phantom used). Based on the empirical formula 
provided in Refs. 27 and 28, a can be calculated to be 0.63 for the estimated boundary 
surrounding. Our derived values of a appear to have a good agreement with this estimated 
actual value. Since a mixed surrounding medium is often used in reality for optical tomography, 
our method offers a practical way to handle the complex experimental boundary conditions. The 
small difference in a values for the single and two target cases are due to the difference in the 
composition of the phantom materials used in each case. With the derived a values, we obtained 
the optimal reconstructed images shown in Fig. 5. As can be seen, both single and two targets 
can be clearly detected quantitatively in terms of the optical property value, size and location of 
the targets. 

In summary, we have developed two schemes that can enhance our established 
reconstruction algorithms. The described normalizing scheme-based algorithm allows for 
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reconstruction of optical property images without the calibration procedure with respect to the 
excitation source. Furthermore, our simulated results show that the normalizing scheme-based 
algorithm is even more resistant to the noise effect than the reconstruction algorithm without a 
normalizing scheme. We have also demonstrated that the use of a simple least-squares 
minimization scheme can provide us the correct boundary conditions coefficient that is important 
for a successful reconstruction. Our simulations and experiments indicate that the concurrent use 
of the two schemes can eliminate the calibration procedure that was previously needed when 
dealing with experimental data. Although we have chosen cw illumination as our test case in this 
study, the proposed schemes can be easily adapted for reconstructions in frequency- or time- 
domain. 

This research was supported in part by grants from the National Institutes of Health 
(NIH) (CA 78334), the Department of Defense (DOD) (BC 980050), and the Greenville Hospital 
System/Clemson University Biomedical Cooperative. 
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Figure Captions 

Figure 1. Images reconstructed by non-normalizing algorithm (Non-NRA) and normalizing 
algorithm (NRA) from simulated data without noise added to the “measured” data, (a) |i a image 

by Non-NRA; (b) |i' image by Non-NRA; (c) p a image by NRA; (d) \i' s image by NRA. 

Figure 2. Images reconstructed by non-normalizing algorithm (Non-NRA) and normalizing 
algorithm (NRA) from simulated data with 5% noise added to the “measured” data, (a) |i a image 
by Non-NRA; (b) |i' image by Non-NRA; (c) p a image by NRA; (d) \l' s image by NRA. 

Figure 3. (a) X 2 as a function of a calculated without noise added to the “measured” data. 
Three sets of optical property values were used for computing X 2 : (1) 0%-off the background 
values (the solid line), (2) 2%-off the background values (the dash-dotted line), and (3) 20%-off 
the background values (the dashed line), (b) The same as (a) except that 5% noise was added to 
the “measured” data. 

Figure 4. X 2 versus a for experimental data when both single and two target configurations 
were used. 

Figure 5. Reconstructed images with experimental data, (a) |i' image for the single target case, 
(b) (l a image for the single target case, (c) |l' image for the two target case, (d) p a image for the 
two target case. 
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Abstract: We have developed a fully automated frequency-domain system for combined optical 
and fluorescence reconstruction imaging. The system is currently being evaluated using tissue 
phantom experiments. 
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OCIS Codes: (170.3880) Medical and biological imaging; (170.5270) Photon density waves; (170.3650) Lifetime-based 
sensing; (110.0110) Imaging systems; (110.6960) Tomography. 

We report the design and development of a near-infrared (NIR) frequency-domain system that is be able to provide 
quantitative optical and fluorescence images of breast tissue. We are presently using phantom experiments to 
evaluate our system. 

Our efforts have been focused on designing a system that can offer high signal-to-noise ratio in order to 
obtain higher measurements accuracy and thus a good contrast of reconstructed optical/fluorescence images. 
Because a multiple source-detector arrangement was used, attention was paid to carefully measure the optical 
transmission over every source-detector path and as well as to the introduction of these corrections in the 
reconstruction algorithm. 

Our imaging system is schematically shown Fig. 1. 



This system employs radio-frequency intensity modulated near-infrared beam from a 60 mW diode laser at 
789nm (Mitsubishi). The bias current is modulated at 100 MHz using an IFR 2032A RF Generator (IFR Americas). 
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The laser beam is sent to the phantom by 16 fiber optic bundles (Fiber Optics, Inc.) with 3 mm inner diameter each, 
coupled with a high precision Melles Griot moving stage. The diffused radiation is received by the second 16 
channel fiber optic bundles and delivered to a PMT (R928 -Hamamatsu Corporation). A second PMT is used to 
record the reference signal. These photodetectors are supplied at a radio frequency modulated current with 0.1-lKHz 
shift The intermediary frequency signal obtained from the PMTs is processed using a DAQ 6032 National 
Instruments board (16 bits, lOOKHz sampling rate). All elements are computer controlled by a GPIB bus. To 
increase the dynamic range of intensities that the PMT can detect, an automated filter wheel with pre-calibrated 
neutral density filters is added to the system. Two low pass filters with a cut off of 5 kHz are used to reduce the high 
frequency noise. AC, DC, and phase shift between reference and probe signals are obtained through the LAB VIEW 
routines. 

For every position of the source, 16 measurements for each detector are made, taking alternatively 100 ms 
samples for fluorescence and nonfluorescence signals. Fluorescence signals are measured using an. 830nm 
interferential filter placed in front of photomultiplier. 

Almost simultaneously (10" 5 sec delay), reference signal is taken in order to measure the phase shift. By 
using a sample and hold scheme, this time can be reduced at 5 ns and so the errors in the measuring of the phase 
shift can be further reduced. FFT is used to obtain real and imaginary parts of the reference and probe signals. The 
total data acquisition time for 256 measurements is about 8 minutes at present time. This acquisition time is long but 
can be reduced by using a higher speed moving stage or by using a PMT array. 

The performance of the system is strongly dependent upon the quality of optical alignment of FO bundles, 
the characteristics of the modulated light and quality of detectors. Our efforts were directed at obtaining a high 
signal to noise ratio (SNR) for every position of the detector. In particular, a high DC intensity and a high 
modulation ratio are very important in order to keep a high SNR. Some optical collimating systems have been used 
to collect all radiation passing through the ND filters and focusing onto the PMT cathodes. 

We have performed measurements for a 50 mm diameter cylindrical solid phantom (0.5% Intralipid +Agar) 
containing a 7.5 mm diameter target cylinder (0.5% Intralipid + ICG dye + Agar). The measured relative AC and 
phase shift for both excitation and fluorescence signals are presented in Fig.2. The intensity and phase shift of 
fluorescence signal are governed by the local fluorescent properties (quantum yield and lifetime). In addition, the 
fluorescent light that is generated from within the tissue and that propagates to the tissue surface suffers yet another 
round of attenuation and phase delay as a result of the optical properties of the intervening tissue. 



Fig. 2 AC rd and phase shift at different relative position of detector. The fluorescent dye used was ICG. 

These measurements have been repeated several times in order to examine the repeatability of the system 
and to calibrate the system. After the system was calibrated, several sets of 256 measurements were obtained for a 
cylindrical configuration. Attention was paid to the alignment of axes of the optical fibers that scan the 16 emitting 
fibers and the 16 receiving fibers. The translation stages that are used to multiplex the fiber optics have a 
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repeatability of approximately 1pm, allowing a very low variation in the measured signals. The repeatability of the 
system is approximately 0.5 degree in phase and 0.2 % in AC and DC signals. 

The measured excitation and fluorescent data, are currently being used to reconstruct both optical and 
fluorescence images. Our reconstruction algorithms have been described previously in Refs. 1 and 2. 


This work was supported in part by the NIH grant (CA 78334) and the DOD grant (BC 980050). 
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Abstract 


We experimentally demonstrate what is believed to be the first simultaneous 
reconstruction of fluorescence lifetime and yield distributions in heterogeneous turbid media 
using absolute frequency-domain measurements. The fluorescence images are obtained by use of 
a finite element based reconstruction algorithm. Experiments are performed using indocyanine 
green (ICG) dye and tissue-like phantoms in both single- and multi-target configurations. 

OCIS Codes: 170.3010,170.3650,170.3830, 170.6280,170.6960 
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Recently the idea of reconstruction-based fluorescence diffusion tomography has been 
proposed and developed for clinical applications such as breast cancer detection and tissue 
functional mapping. 1 ' 5 This new imaging approach relies on the fact that the lifetime and yield of 
fluorophores in tissue can potentially provide tissue functional information such as tissue 
oxygenation, pH, and glucose. 6 It is also based on the fact that the fluorophore will preferentially 
accumulate in tumors, hence providing the best sensitivity for cancer detection. 7 In this type of 
indirect imaging method, an effective reconstruction algorithm is crucial which allows for the 
formation of a spatial map of the lifetime and yield of heterogeneous fluorophore distributions. 
To date fluorescence image reconstructions are largely limited to simulated data, based on linear 
and nonlinear algorithms. 1 * 5 Chang et al. 2 attempted to obtain a qualitatively good image of the 
fluorophore concentration only from dc measurements. In this Letter we present successful 
quantitative simultaneous reconstruction of both lifetime and yield images from frequency- 
domain measurements using indocyanine green (ICG) dye and tissue-like phantoms. To our 
knowledge, the results shown here are the first experimentally reconstructed images for which 
simultaneous recovery of both lifetime and yield profiles has been achieved with absolute ac 
excitation and fluorescence data. We show the ability of our algorithm to reconstruct images in 
which single- and multi-target are embedded in a circular background. 

We have used in this study the finite-element based algorithm for reconstruction that has 
been described in detail in Ref. 5. The algorithm uses a regularized Newton method to update an 
initial (guess) fluorescent property distribution iteratively in order to minimize an object function 
composed of a weighted sum of the squared difference between computed and measured data. 
Here we present a brief overview for context. 
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In frequency-domain, it is known that propagation of both excitation and fluorescent 
emission light in tissues or turbid media can be described by the following coupled diffusion 
equations: 2 ' 5 


V.[D,(r)V<D,(r,m)] 
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M r >- 

(A,_ <r) — 


ICO 

c 

ico 

c 


<J> x (r,co) = -S(r,co) 


( 1 ) 

( 2 ) 


where <J> xm is the photon density for excitation (subscript x) or fluorescent light (subscript m), 
D x m is the diffusion coefficient, p. a n is the absorption coefficient due to contributions from 
both non-fluorescing chromophores and fluorescent dye, |i a is the absorption coefficient for 

the excitation light due to contribution from fluorescent dye, © is the modulation frequency, c is 
the velocity of light in the medium, andi] and x are the fluorescent quantum yield and lifetime, 
respectively. S(r,co) is the excitation source term in (1). Note that a single-exponential 
fluorescence decay has been assumed in the source term for fluorescent light (right-hand-side of 

(2)); multi-exponential time decay can be handled by a simple extension. The diffusion 
coefficient can be written as D xm = 1/3 (|_l^^ +(l a a ) where ^ is the reduced scattering 

coefficient. In this study a point excitation source, S = S 0 8(r-r 0 ), and Type HI boundary 

conditions (BCs): -D xm V<l) xm - h = ad> xm , are used, where S 0 is the source strength, n is the 

unit normal vector for the boundary surface and a is a coefficient that is related to the internal 
reflection at the boundary. Both S 0 and a are determined directly from the heterogeneous 

measurements by use of a data pre-processing scheme that has been reported elsewhere. 8 

Making use of finite element discretization, we can obtain two matrix equations for Eqs. 
(l)-(2) and realize other derived matrix relations through differentiation, which yield a set of 
equations capable of inverse problem solution: 
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where the elements of matrices [A x ] and [A m ] are respectively 

( a x, m )ij =<- D x m V\jrj -V\|/, -(|i a n and the entries in column vectors {b x , ra } and 
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^ ^x.m = {(^x.tJi’^x.iJi,' where ( ) indicates integration over the problem 

domain, and <J> x m , x and ti have been expanded as the sum of coefficients multiplied by a set of 
locally spatially-varying Lagrangian basis functions \|/j, \|/j and \j/ k . j expresses integration 
over the boundary surface where Type IE BCs have been applied.. (O xm )j is the photon density 
at node i, N is the node number of a finite element mesh and M is the boundary node number. % 
expresses x or T|. and Ojf are the observed and computed fluorescent data, respectively. In 
fluorescence lifetime tomography, the goal is to update the x and t| distributions through the 
solution of Eqs. (3)-(6) so that a weighted sum of the squared difference between computed and 
measured data can be minimized. 

Our experimental setup used was an automated multi-channel frequency-domain system 
that was described in detail elsewhere. 9 The system employed a radio-frequency intensity- 
modulated near-infrared beam. The laser beam at 785 nm was sent to the phantom by 16 fiber 
optic bundles coupled with a high precision moving stage. The diffused radiation was received 
by another 16 channel fiber optic bundles and delivered to a PMT. A second PMT was used to 
record the reference signal. These PMTs were supplied at a radio-frequency modulated current 
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with 0.1-lKHz shift. The intermediary frequency signal obtained from the PMTs was processed 
using a National Instruments board. To increase the dynamic range of intensities that the PMT 
can detect, an automated filter wheel with pre-calibrated neutral density filters was added to the 
system. For every source position, 16 measurements for each detector were made, taking 
alternatively 100 ms samples for sample and reference signals. Fluorescence signals were 
obtained through an 830nm interference filter placed in front of the detection PMT. ac intensity 
and phase shift between reference and sample signals were obtained using FFT Labview 
routines. The total data collection time for 256 measurements was 8 minutes. 

In our experiments we used a 50 mm diameter solid phantom (0.5% Intralipid +India 
ink+Agar) as the background medium, making the absorption coefficient (i a = 0.0005 / mm, and 

the reduced scattering coefficient p' = 0.5/mm. Both single- and two-target configurations were 

used. The solid targets contained ICG dye with different concentrations (0.5 or 1 pM) plus 0.5% 
Intralipid and Agar (perfect uptake of the dye was considered here, i.e., no dye was present in the 
background). The single target had 12.6 or 7.8 mm in diameter, while both targets had 10.3 mm 
for the two-target case. The 2D finite-element mesh used had 249 nodes and 448 elements for 
both forward and inverse solutions. 

Figs. 1-4 show the reconstructed images of both t and rj and their property profiles 
along a cut-line crossing the background and target regions for one off-centered target having 
different target size (Figs. 1 and 2) and different dye concentration (Figs. 1 and 3) and two off- 
centered targets (Fig. 4). As can be seen, all the targets can be clearly imaged. We also noticed 
that the recovered images not only are quantitative with respect to the location and size of the 
target but also quantitatively resolve the values of the lifetime and yield distributions. From Figs. 
1-4, we can see that the sharp edges of the target for T| images can be recovered, while the target 
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edges are generally broadened for x images. However, the values of x images are much more 
accurately reconstructed than that of r\ images, as evidenced by the property profiles shown at 
the bottom of each figure. This is not surprising since normally one would not be able to 
reconstruct a pure t| image from fluorescence diffusion tomography. Instead, one would obtain 
an image of the product of the yield and absorption coefficient due to the dye. 1-5 Here we have 
assumed this absorption coefficient to be known as a constant for the entire problem domain. 
This approximation could be the primary contributor to the inaccuracy of tj image recovery. It is 
also interesting to note that a finite value of x (larger than the x value in the target) for the 
background region is always reconstructed in all the cases studied, although there is no 
fluorophore distribution in the background. 

In summary, we have presented simultaneous reconstruction of both lifetime and yield 
images in turbid media using ac excitation and fluorescent data. While an ideal, perfect uptake of 
fluorescent dye in the target has been assumed, this study has clearly demonstrated the feasibility 
of fluorescence lifetime imaging in turbid media using the reconstruction approach described 
here. To add fluorophores into the background would primarily challenge our hardware system 
since the fluorescent signals would be weakened in this case. However, we do not anticipate 
considerable difficulty to obtain enough signal-to-noise ratio for reconstruction when the 
fluorophores are present in the background since the dye uptake ratio between the tumor and 
normal tissue can be as high as 10:1. 4,5,7 In this regard, we plan to conduct experiments in the 
near future. 

This research was supported in part by grants from the National Institutes of Health 
(NIH) (CA 78334), the Department of Defense (DOD) (BC 980050), and the Greenville Hospital 
System/Clemson University Biomedical Cooperative. 
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Figure Captions 

Figure 1. Reconstructed x ((a), top) and r\ ((b), top) images and their property profiles along a 
cut-line (bottom) of a circular target embedded in. the background. The target contained 1 |iM of 
ICG dye and its size was 12.6 mm in diameter. No dye was included in the background. Note 
that the exact values of x and r| for ICG are 0.58 ns and 0.016 in water, respectively 7 ; the exact 
values of x and r| for the background were assumed to be 2.0 ns and 0.0, respectively. 

Figure 2. Reconstructed x ((a), top) and r| ((b), top) images and their property profiles along a 
cut-line (bottom) of a circular target embedded in the background. All the parameters used here 
were the same as in Fig. 1 except that the target had 7.8 mm in diameter. 

Figure 3. Reconstructed x ((a), top) and r\ ((b), top) images and their property profiles along a 
cut-line (bottom) of a circular target embedded in the background. All the parameters used here 
were the same as in Fig. 1 except that the target contained 0.5 |iM of ICG dye. 

Figure 4. Reconstructed x ((a), top) and r\ ((b), top) images and their property profiles along a 
cut-line (bottom) of two circular targets embedded in the background. All the parameters used 
here were the same as in Fig. 1 except that the targets had 10.3 mm in diameter. 
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